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Using three-dimensional numerical simulations, we demonstrate the growth saturation of an un- 
stable thin liquid film on micropatterned hydrophilic-hydrophobic substrates. We consider different 
transverse-striped micropatterns, characterized by the total fraction of hydrophilic coverage and 
the width of the hydrophilic stripes. We compare the growth of the film on the micropatterns to 
the steady states observed on homogeneous substrates, which correspond to a saturated sawtooth 
and growing finger configurations for hydrophilic and hydrophobic substrates, respectively. The 
proposed micropatterns trigger an alternating fingering-spreading dynamics of the film, which leads 
to a complete suppression of the contact line growth above a critical fraction of hydrophilic stripes. 
Furthermore, we find that increasing the width of the hydrophilic stripes slows down the advancing 
front, giving smaller critical fractions the wider the hydrophilic stripes are. Using analytical approx- 
imations, we quantitatively predict the growth rate of the contact line as a function of the covering 
fraction, and predict the threshold fraction for saturation as a function of the stripe width. 



I. INTRODUCTION 

The growth of patterns in forced thin liquid films on 
solid substrates is an ubiquitous process of both bio- 
logical and technological relevance which has attracted 
the attention of the scientific community over the last 
decades [l|, Q . The manipulation of thin films is impor- 
tant in microfluidics, where the increasing miniaturiza- 
tion of fluidic devices demands the development of new 
and simple ways to handle small volumes of liquid. Due 
to their smaller internal hydrodynamic resistance and 
increased surface area, thin films can be used to over- 
come problems associated to other geometries, such as 
microchannels Q. Both experimental and theoretical 
studies indicate the possibility of using clever structur- 
ingof microfluidic devices to orient the growth of films 
J34al or to control the motion of interfaces in confined 
microfluidic chambers (lOj . 

Still, a remaining challenge is how to control the intrin- 
sic growth of thin liquid films in the microscale. When 
a thin film of liquid is forced to spread on a dry solid 
substrate, there is a net accumulation of fluid close to 
the contact line, which corresponds to the apparent posi- 
tion where the liquid, the solid, and the surrounding gas 
meet. As a result, the film profile bulges near its leading 
edge, generating a capillary ridge. For films driven by 
a body force, such as a constant pressure gradient, vari- 
ations in the size of the ridge along the leading edge of 
the film have the effect of increasing the hydrostatic pres- 
sure locally. This introduces a driving mechanism that 
amplifies small perturbations to the contact line, which 
are relaxed by the surface tension [ll[ . Such competition 
of drivin g an d restoring forces gives rise to a linear in- 
stability [13 - [l4| . which, far from a flat advancing front, 



leads to the formation of non-linear structures, such as 
the familiar water rivulets observed in the shower. 

Given the strong interaction with the solid, it is un- 
surprising that the instability is affected by the wetting 
properties of the liquid [15l-[l8l] . Wetting interactions de- 
termine the equilibrium contact angle of the fluid-fluid 
interface with the solid. As a consequence, the size and 
shape of the capillary ridge, which control the driving 
destabilizing force, depend on wetting. For hydrophilic 
substrates, the contact angle of the interface profile is 
small, making the capillary ridge relatively thin. The 
instability is thus weakened Q. Remarkably, if the size 
of the ridge is sufficiently small, the restoring capillary 
pressure is strong enough to balance the driving force af- 
ter the early stages of growth. The interface thus stops 
growing, but still retains a curved saturated shape rem- 
iniscent of a sawtooth pattern (l6j . In contrast, for hy- 
drophobic substrates the net accumulation of mass at the 
ridge is large enough to outweigh the effect of surface ten- 
sion, and the emerging film patterns grow steadily, with 
a shape that resembles that of fingers [17l[l8|. Crucially, 
the lateral lengthscale of the patterns, A max , is of the 
order of the thickness of the thin film. Therefore, in a 
microfluidic device one expects that the lengthscale of 
emerging patterns is also of the order of microns. 

Our aim in this paper is to demonstrate the exploita- 
tion substrate heterogeneity to control the growth of the 
contact line, up to complete growth saturation. This is 
appealing to many situations in which fingering is unde- 
sirable, for example in coating processes and in 'chemi- 
cal channeling' J3(, where the desired lateral lengthscale 
of the film needs not to be constrained by the intrin- 
sic lengthscale of the instability. We shall thus analyze 
the effect of a simple, yet effective, substrate pattern 
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that consists of hydrophilic-hydrophobic stripes oriented 
transversely to the direction of flow, as shown in Fig. [T] 
As we shall see below, the lengthscale of the micropat- 
tern is important, and has a significant impact already at 
scales comparable to the micron-sized films. Still, and de- 
spite its simplicity, such apattern has not yet been stud- 
ied. In a previous study [9j|, we have demonstrated that 
the growth of the thin film can be tuned by patterning the 
solid substrate with transverse or checkerboard arrange- 
ments of hydrophilic-hydrophobic domains. Previous ex- 
perimental and numerical studies of flow patterning [3-0] 
have mainly focused on the effect longitudinal tracks of 
varying flow resistance Q or wetting angle (H?}- The 
overall experiments and numerical simulations evidence 
is that substrate heterogeneity can be used to handle thin 
films, through preferential orientation over hydrophilic 
or low-resistance domains in the substrate, or to par- 
tially control its growth using checkerboard hydrophilic- 
hydrophobic patterns. As we shall show in the following, 
the proposed transverse pattern introduces new physical 
mechanisms, particularly the transverse spreading of the 
film, which lead to the desired saturation of the contact 
line. 

Our approach to this problem will be to perform 
Lattice-Boltzmann (LB) simulations [l9[ of Navier- 
Stokes hydrodynamics coupled to a phase-field binary 
fluid model [20j to account for the interface dynam- 
ics. Previous theoretical approaches to model forced thin 
films on chemically heterogeneous substrates have mainly 
consisted of sharp-interface formulations based on the 
Stokes equations Most of these studies use the 

so-called thin-film equations, which are appropriate to 
model the fluid dynamics in the limit of small interface 
slopes. Within this framework, a common procedure to 
regularize the viscous dissipation singularity is to replace 
the contact line by a thin precursor film that extends be- 
yond the leading front. The effect of the heterogeneity 
can then be included by choosing a particular model for 
the precursor. For example, an heterogeneous disjoining 
pressure model that couples to the precursor film Q has 
been used to study the motion of the film on longitudinal- 
striped and disordered substrates, as well as on ordered 
and disordered isolated spots Q. Alternatively, a sim- 
ilar model consisting of a spatially-varying thickness of 
the precursor film to study the motion of the film over 
chemical patterns has also been used 0, Q . 

Our aim in this work is to study the effect of a sharp 
contrast between hydrophilic and hydrophobic domains 
in the proposed micropattcrns. This corresponds to sit- 
uations in which the static and dynamic contact angles 
need not to be small. Such regime falls beyond the limit 
of applicability of the thin-film equations, demanding the 
resolution of the full three-dimensional dynamics. This 
is a complicated problem from the classical point of view, 
as it involves the motion of two free boundaries: the film 
surface and the contact line. Within sharp-interface for- 
mulations, boundary integral methods could be used to 
address this problem, as has been done to model the effect 



of chemical heterogeneities on two-dimensional fluid nan- 
odrops plj . However, the mesoscopic approach that we 
take here constitutes a reliable three-dimensional Navier- 
Stokes solver that allows for both small and large contact 
angles, necessary to model hydrophilic and hydrophobic 
interactions. Because of the diffuse approximation to de- 
scribe the fluid interface, the method circumvents the 
free-boundary problem, regularizes the viscous dissipa- 
tion singularity at the contact line, and deals naturally 
with merging and pinch-off events, which otherwise need 
a specific rule in sharp interface approximations. We 
have carried out a validation of our model in a previous 
study, for the case of driven thin films on homogeneous 
solid substrates [22T |. 

The rest of this work is organized as follows. In 
Section HQ we introduce the mesoscopic model and the 
Lattice-Boltzmann algorithm to perform the numerical 
simulations. Wc will complement our numerical results 
with kinematical calculations for the propagation of the 
interface. As wc shall see from our numerical results, 
presented in Section [TV] the main effect of the transverse 
hydrophilic stripes on the forced film is to favor its lateral 
spreading on hydrophilic domains, thus slowing down the 
leading edge of the film. The same lateral spreading ul- 
timately leads to a speed-up of the trailing edge of the 
film. The net result is that the growth rate of the film 
decreases with increasing fraction of hydrophilic stripes 
in the substrate, until the film saturates above a critical 
value of this fraction. As we detail in the discussion and 
conclusions of this work, presented in Section [Vj we pre- 
dict that a relatively small fraction of hydrophilic stripes 
is sufficient to suppress the growth of the contact line in 
experimental realizations. Our results thus constitute a 
promising step towards the control of interface growth in 
microfluidic systems. 



II. MODEL 

In this paper we shall perform numerical simulation of 
Navier-Stokes level hydrodynamics coupled to a phase- 
field representation of the viscous fluid film and surround- 
ing low-viscosity phase. Such modeling is very useful to 
address the dynamics of immiscible fluids, as it has the 
advantage of circumventing the free-boundary problem 
associated to sharp-interface representations, as for ex- 
ample in boundary integral methods. Furthermore, the 
phase-field model regularizes the viscous dissipation sin- 
gularity at the contact line through a diffusive mecha- 
nism acting at small scales (of the order of the size of 
the interface). One thus does not need to specify addi- 
tional boundary conditions at the contact line. In order 
to integrate the model equations, we will use a Lattice- 
Boltzmann algorithm. 
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FIG. 1: Schematic representation of the system. A thin film of 
initial thickness h is forced along the x direction by the action 
of a constant body force, /, on a solid substrate composed of 
hydrophilic and hydrophobic transverse stripes. The scale of 
the thin film thickness and of the width of the stripes is of 
the order of microns. 



A. Hydrodynamics and Phase-Field Model 

A schematic depiction of the forced thin-film geometry 
is presented in Fig. [TJ We consider two immiscible New- 
tonian fluids of viscosities rji and 772 and densities pi and 
P2, respectively. The thin film, corresponding to fluid 1, 
is forced on a solid substrate composed of hydrophilic and 
hydrophobic regions. In order to distinguish between the 
two liquids, we introduce a phase field characterized by a 
concentration variable, or order parameter, <f>. The local 
value of the order parameter thus describes each phase, 
and varies between two volume values across a diffuse 
intcrfacial region. 

In equilibrium, the free energy of the system can be 
written as a functional of the order parameter field, 0(r) 
and the local density field, p{r), 



HP,<P} = J Av{u(^p) + ^{V^f)+ J dSfs^s). 

(1) 

The first integral in Eq. (JXJ) is composed by the vol- 
ume contributions to the free energy. These consist of a 
double-well potential with an ideal gas term, U(4>, p) = 
A<f> 2 /2 + Bcf) 4 /4+ (l/3)plnp, and a square-gradient term 
that penalizes variations of the order parameter, thus 
allowing for the formation of a diffuse interface. The pa- 
rameters A, B and k control the equilibrium values of 
the order parameter ±0 cq 
interface, £ 

sion, 7 = y/— 8kA a /9B 2 . The second integral in the free 
energy corresponds to the contribution of fluid-solid in- 
teractions. Here we chose a Chan model, fs{4>s) = C0s, 
which assigns the free surface energy fs according to the 
local value of the order parameter at the boundary, <j)g. 
The parameter C controls the equilibrium contact angle, 
6L, of the fluid-fluid interface in contact with a flat solid 



± \J —A/ B, the size of the 
\/—k/2A and the fluid-fluid surface ten- 
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+ [1 + C(-kAY 



(2) 

The order parameter and density fields contribute to 
the chemical potential, p, and pressure tensor, P a /3, of the 
system. These are obtained by taking functional deriva- 
tives of the free energy [23| and have the following ex- 
pressions, 



p = A<t> + B<f) 3 - kV 2 



(3) 



P 



a/3 



\A<f?- 
-n(d a (j))(d cj)). 



(4) 

Out of equilibrium, the dynamics of the density, p, 
velocity, v, and order parameter, </>, fields are given 
by the continuity and Navicr-Stokcs equations, plus a 
convection-diffusion equation for the phase field, i.e., 



| + v.(^ = o, 



(5) 



dv 



+ (v- V)v = -VP - <t>Vp + r)V 2 v + f, (6) 



and 



30 
dt 



V<p = MV 2 p, 



(7) 



respectively. 

The 0— dependent term in the Navier-Stokes equations 
arises from the concentration dependence of the pressure 
tensor in Eq. (|4]), giving a 'chemical' force density, —(j)Vp, 
which plays a similar role as the pressure gradient term 
in the volume of each phase. The remaining terms in 
the Navier-Stokes equations are the viscous friction term, 
whose strength is controlled by the local viscosity, 77, and 
to the body force term, /. 

The evolution of the order parameter is dictated by 
Eq. (0, which contains an advective term caused by 
the underlying velocity field and a diffusive term, whose 
strength is controlled by the mobility parameter M . 



B. System Geometry and Boundary Conditions 

In order to ensure the formation of steady patterns, 
we choose a constant flux configuration for the geome- 
try of the system. We set a rectangular simulation do- 
main of linear dimensions L x W x H, in the x, y and 
z directions, respectively. The solid-fluid interface, S, 
is parallel to the x — y plane and located at z = as 
shown in Fig. [TJ The thin film is initially at rest and 
occupies the volume Vq = Lq x W x h c . We fix the 
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force term as / = ^f x (<fi + l)x, so fluid 1 is forced uni- 
formly while fluid 2 is left to evolve passively. This sets 
up a flow of mean velocity u = h^.f x /3rji within the film. 
At x = 0, a constant flux in the x direction is thus en- 
sured as long as there is a reservoir of fluid available 
to compensate for the downstream motion of the film. 
This is achieved by imposing d x p(x = {0,L},?/, z) = 0, 
d x v{x = {0, L},y, z) = and d x (j)(x = {0,L},y 7 z) = 0. 
Using this approach, we do not find any appreciable vari- 
ations of the mean velocity of the film or the film thick- 
ness close to the boundary. 

We fix periodic boundary conditions in the y direction, 
i.e., p(x,y = W,z) = p(x,y = 0,z), v(x,y = W,z) = 
v(x,y = 0,z), (j>(x,y = W,z) = (j>{x,y = 0,z). At 
the upper boundary, a shear-free boundary condition is 
imposed by fixing d z v{x, y, z = H) = 0, while vanish- 
ing density and concentration gradients normal to the 
boundary are fixed by imposing d z p(x,y,z = H) = 
and d z 4>{x, y,z = H) = 0. 

The velocity profile is subject to stick boundary and 
impenetrability conditions at the solid- fluid surface, i.e., 

v\s = 0. (8) 



Eq. (|T0|) corresponds to the Young-Laplace condition 
wit h R being the local radius of curvature of the interface. 
Thus, the concentration dependent term in the Navier- 
Stokes equation eliminates the need of using Eq. (fTU)) 
explicitly. Accordingly, the continuity of the tangential 
stress is recovered in the sharp interface limit. 

In order to present our results in a natural scale, we 
choose units according to the classical thin-film theory in 
the lubrication limit p], [|[ . For microscales inertial effects 
are negligible, and the relevant control parameter is the 
capillary number, Ca = 77117/7, which measures the ratio 
of viscous to capillary forces. A suitable unit system is 

(*>' = —, and * (II) 
u t c 

where x c = h c (3Ca)~~s , and t c — x c /u. 



In macroscopic representations, imposing the stick 
boundary condition for a moving interface in contact with 
a solid causes a spurious divergence of the viscous stress 
as one approaches the contact line (24[. This paradox 
is readily avoided by the phase-field model, through the 
diffusive term in Eq. (J7J, which allows for a slip veloc- 
ity that increases as one approaches the contact line [25[ , 
even when Eq. ([8]) is enforced. Diffusion is driven by 
variations of the chemical potential profile close to the 
contact line, over a typical lengthscale Id, which scales as 

An evolution equation for the order parameter should 
be specified at the solid-fluid interface (see [25| for a de- 
tailed derivation). In this paper we shall impose the nat- 
ural boundary condition 



n ■ V</>|< 



nd(/)s' 



(9) 



This gives an instantaneous relaxation of the order pa- 
rameter gradient to its surface value at the wall, which 
corresponds to quasi-equilibrium dynamics for the inter- 
face slope at lengthscales comparable to the size of the 
interface (23|- By choosing the wetting parameter C in 
Eq. ([2]), we specify the local value of the equilibrium con- 
tact angle and hence the desired hydrophilic-hydrophobic 
pattern. 

The phi— dependent term in Eq. ([6]) is most important 
at the fluid-fluid interface, where it plays the role of a 
Laplace pressure. This behavior can be verified by in- 
tegration of Eq. over the length of the interface (c.f. 
|28( ) whence one recovers a jump in the normal stress 



R' 



(10) 



C. Parameter values 

The model presented in the previous section is specified 
by a set of parameters, pi, p2, J?i, T}%, 4>cq, 7, £ and M. 
Additionally, the fluid wetting properties are specified by 
the static contact angle, 9 e , which depends on C . In the 
following, we specify the values of these parameters used 
in our simulations (all in lattice units). 

The equilibrium values of the order parameter are fixed 
to 0* = ±1 by choosing A = — B. To specify the viscosity 
of each fluid we fix the mean viscosity (77) = 0.1 and the 
viscosity contrast, 5-q = 0.9. The local viscosity, r)(r), 
is then set as r/i = (77) (1 + 5rf) if <f> > (viscous film) 
and 772 = + Srj) if <f> < (surrounding low- viscosity 
phase) . The density is set to the same mean value in each 
phase, pi = p2 = 1. The fluid- fluid surface tension is 
fixed to 7 = 2.3 x 10~ 3 . The interface width is set to £ = 
0.57, which corresponds to an effective transition zone 
of ~ 5 lattice spacings for the order parameter profile 
between the volume values in equilibrium. 

In order to ensure that the contact line propagates with 
a velocity comparable to the velocity of the film [22j , we 
fix the mobility to M = 10.0. 

The equilibrium contact angles used throughout this 
study correspond to 9 e = 0° for the hydrophilic substrate 
and 8 e = 90° for the hydrophobic substrate, which corre- 
sponds to C = and C = 1.7 x 10~ 3 , respectively. 

The capillary number is set to Ca = 0.41, which follows 
from choosing the forcing f x to give a mean velocity value 
for the film, u = 0.005. This value ensures the stability 
of the LB scheme and a small compressibility (density 
variations of the order of ~ 



III. LATTICE BOLTZMANN METHOD 

We use the Lattice-Boltzmann (LB) implementation 
introduced in detail in Ref. 0, [29| in the context of spin- 
odal decomposition of binary fluid mixtures, and which 
we have used in a previous work where we have fo- 
cused on the dynamics of unstable thin films in homo- 
geneous solid substrates. In LB, the dynamics are intro- 
duced by two sets of velocity distribution functions, fi 
and gi, which evolve in time according to the discrctized 
Boltzmann equations, 




/<(?+S,t + l)- fi(r,t) 



(fi-m+Fi, (12) 



and 



g l (r + Ci,t+1)- gi(f, t) = (g t 

Tn 



(13) 



Space is discretized as a cubic lattice where nodes are 
joined by velocity vectors, c). Here we use the D3Q19 
model, which consists of a set of 19 velocity vectors in 
three dimensions [23|]. Hence, fi and gi are proportional 
to the number of particles moving in the direction of c\. 
Eqs. (fT2j) and (fT3|) are composed of two steps. First, 
the distribution functions arc relaxed to their equilib- 
rium values, represented by f? q and g eq , with relaxation 
timescales tj and r g . The term F/ is related to the ex- 
ternal forcing. Following this collision stage, distribution 
functions arc propagated to neighboring sites. 

The mapping between LB scheme and the hydrody- 
namic phase-field model is done by defining the hydro- 
dynamic variables through moments of the fi and gi. 
The local density and order parameter are given by 
J2i f i = P an d J2i 9i = <Pt while the fluid momentum 
and order parameter current arc defined as fi&i = pv 
and J2i 9iCi = 4>v. 

Local conservation of mass and momentum is enforced 
through the conditions J2i ft = J2i /f 9 = P> I2i 9i = 
YLidT = 0, EiM = Ei/r^ = Pv and Y<i9& = 
J2i 9\ q Ci = 4>v. In equilibrium, the pressure tensor and 

chemical potential are defined as f^ q Ci&i = pvv + Pt 

and 9i 9 <k<k = M/j,S + cfrvv. 

To recover the equilibrium behavior of the phase-field 
model, the equilibrium distribution functions and the 
forcing term are expanded in powers of €?[30|| . i.e., 



ft, 



At 



+ Sv ■ Ci + -vv 



dCi - -v 



= f 

G : dc, 



g eq = puj u ( A q + 3v ■ Q + ^-vv : SiCi - ^v 2 + G 9 : c,c 



and 



F =4uu\l 



2t, 



f ■ ci(l + v-Ci) -v- f 



FIG. 2: Geometry of the transverse stripe pattern. The pat- 
tern is composed of alternated hydrophilic (grey) and hy- 
drophobic (white) stripes of lengths lu and 1$, respectively. 
The panel shows a schematic representation of the contact 
line (solid line), which spreads on the hydrophilic stripes up 
to a width w^W before advancing on the next hydrophobic 
stripe. Spreading causes that a series of dry spots are left on 
the hydrophobic stripes, which are eventually covered by the 
trailing edge of the film. The dotted lines depict the thin film 
profile, which has a typical thickness ft? on a wet hydrophilic 
stripe, and a maximum thickness h t at the trailing edge on a 
hydrophobic stripe. At the grey zones 9 e = 0° while at the 
white zones 9 e — 90°. 



Here, v stands for the three possible magnitudes of the 
Ci set. The values of the coefficents are then fixed by 
the definition of the moments of f^ q and g eq ', and by the 
symmetry of the lattice. Coefficient values for the D3Q19 
LB scheme are ojq = 2/9, w\ = 1/9 and uj^ = 1/72; 

A f = 9/2 - 7/2TrF, A[ = A f ^ = 1/pTvP and & = 

9/(2p)P - 3<5TrP; A 3 = 9/2 - 21/2M/X, A? = A 9 ^ = 

3Mp/p and G 9 = 9/(2p)Mp(I - I), where I is the unit 
matrix. 

The hydrodynamic phase-field model, given by 
Eqs. ([5]), ((6j) and (J7J, up to second order corrections in 
the velocity, can be recovered by performing a Chapman- 
Enskog expansion of Eqs. (JT2} and (|T3]> (30|. The LB 
scheme maps to the hydrodynamic model through the 
relaxation timescales, i.e., rj = (2tj — l)/6 and M = 
(r g — 1/2)M, and through the body force / = pg. 

Solid boundaries in the Lattice-Boltzmann method are 
implemented by means of the well known bounce-back 
rules [30l l3ll ]. In the lattice nodes that touch the solid, 
the propagation scheme is modified so the distribution 
functions are reflected to the fluid rather than absorbed 
by the solid. As a consequence, a stick condition for 
the velocity is recovered approximately halfway from the 
fluid node to the solid node. 



IV. RESULTS 

We want to assess if a small fraction of hydrophilic do- 
mains can induce a significant reduction of the growth 
of the contact line. We thus characterize the trans- 
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FIG. 3: Contact line evolution on transverse stripe patterns 
for hydrophilic stripes of length lh/xc = 1.3 at different cover- 
ing fraction, /. The first panel in each figure shows the initial 
perturbation at t/t c = 0, while the subsequent two panels 
show the contact line at time intervals of At/t c — 67. For 
comparison, the last panel in (a) and (b) shows the finger 
that grows on a homogeneous hydrophobic substrate (/=0) 
at t = 2At/t c . For both values of /, the contact line grows as 
a finger on the hydrophobic stripes, and spreads transversely 
on the hydrophilic domains, (a) For a small fraction of hy- 
drophilic stripes, / = 0.11, the contact line grows at a smaller 
rate than the one measured on a homogeneous substrate, (b) 
Increasing the fraction to / = 0.55 has the effect of decreas- 
ing the growth rate dramatically, making the length of the 
contact line substantially smaller than the length of the un- 
perturbed finger. As seen by comparing the position of the 
contact line in third panel of (a) and (b), the reduction of the 
growth rate arises both from a slow down of the leading edge 
and a speed up of the trailing edge of the contact line. 



verse stripe arrangement by the width of the hydrophilic 
stripes, lh, and the fraction of the substrate that is cov- 
ered by them, / = lh/(lh + l<p)i where 1$ is the length of 
the hydrophobic stripes, as depicted in Fig. [2j 

In a previous study @, we have considered the mo- 
tion of thin films on homogeneous substrates. We have 
found that the front develops as a finger if the substrate 
is hydrophobic (6 e = 90°), with an inter- finger spacing 
given by h max (Q e = 90°) ~ 8x c . In contrast, for a hy- 
drophilic substrate the front develops as a sawtooth with 
a typical lateral lengthscale of A max (8 e = 0°) ~ 14x c . 
Due to the periodicity of the emerging structures, it 
is sufficient to consider simulation domains of width 
W = A max (8 e = 90°) ~ 8x c . Larger system sizes obey- 
ing W = nA max (8 e = 90°), with n being an integer, 
only give rise to additional identical fingers, whose effect 
is already taken into account by the periodic boundary 
conditions in our simulations. For non-integer n we ex- 
pect a non-linear competition of emerging fingers, which 
we have observed previously Q. However, this compe- 



tition should not modify significantly the main effect of 
the transverse-stripe pattern. 

In the following, we will explore the dynamics of the 
front when varying lh and /. We study the effect of the 
occupation fraction in the range < / < 1, while for the 
length of the stripes we choose the values l%jx c = 1.3, 
lh/%c = 4.0, and lh/x c = 6.7. As before, we consider 
the evolution of a single mode in systems whose width 
is chosen to be close to the most unstable wavelength 
predicted for the hydrophobic substrate. 

Figures[3la) and[3^b) show the evolution of the contact 
line for / = 0.11 and / = 0.55 respectively, for a stripe 
width lh/x c = 1.3 at three subsequent times, ta/t c = 0, 
ti/tc = 67 and t%lt c = 134. For comparison, the last 
panel at the bottom of (a) and (b) shows the interface 
profile for a homogeneous hydrophobic substrate (/=0) 
at t-iltc = 134. The main feature in these patterns is that 
the growth rate of the contact line decreases strongly as 
the fraction of hydrophilic stripes is increased. For the 
value of lh considered in Fig. [3l the growth rate actually 
vanishes at a critical fraction f c ~ 0.59. We first notice 
that for both patterns shown in Fig. [3l the leading edge 
of the contact line grows as a finger when it is in con- 
tact with a hydrophobic stripe. The finger grows until 
it touches a hydrophilic stripe, where it starts to spread 
transversely. As spreading takes place, the leading edge 
continues advancing over the hydrophilic domain, until 
it reaches the next hydrophobic stripe, where it devel- 
ops as a finger again. A closer inspection of Fig. [3]Ja) 
shows that a series of spots are left uncovered momen- 
tarily on the hydrophobic stripes after the leading edge 
has passed by. This is a consequence of the transverse 
spreading process that occurs every time the leading edge 
passes over a hydrophilic stripe. The second and third 
panels in Fig. EJa) show that the trailing spot is eventu- 
ally covered by the film, and that only after the trailing 
spot has disappeared, the following spot starts to be cov- 
ered. This mechanism holds for larger /, as can bee seen 
in Fig. OJb). The global result is that the transverse 
pattern induces a slow down of the leading edge, as can 
be seen by comparing the contact line with the finger 
that grows on a homogeneous substrate. This effect is 
stronger for / = 0.55 than for / = 0.11. Additionally, 
the transverse pattern induces a speed up of the trailing 
edge, and effect that can also clearly be seen in Fig. [3] 
Again, this effect is stronger the larger / is. 

The combination of the speed up of the trailing edge 
and the slow down of the leading edge causes the contact 
line growth rate, L, to decrease with increasing fraction of 
hydrophilic stripes. Figure [5] shows a plot of L as a func- 
tion of /. The limits in this plot correspond to the finger 
(/ = 0) and sawtooth (/ = 1) solutions on homogeneous 
substrates. For small /, the growth rate converges to the 
growth rate of the unperturbed finger, while increasing / 
has the effect of decreasing the growth rate as explained 
above. We find that L vanishes at a critical fraction, 
fc(lh), which is considerably smaller than unity. This 
means that the transverse pattern induces the saturation 
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l h = 1.3 
h. = 4.0 
l h = 6.7 
pred. wi, — 0.58 
pred. Wh = 0.60 
pred. = 0.64 



FIG. 4: Global contact line growth rate as a function of the 
fraction of hydrophilic stripes for different hydrophilic stripe 
lengths on transverse striped patterns. The growth rate L, 
decreases as the fraction of hydrophilic stripes, /, increases. 
For fixed /, increasing the length of the hydrophilic stripes, 
lh, has the effect of decreasing the growth rate. Continuous 
lines correspond to the kinematical theoretical prediction, for 
which the dependence on lh is contained in the parameter Wh, 
which we measure from simulations. The plot shows that the 
growth rate vanishes at a fraction that is smaller than unity. 



of the contact line by caused by transverse spreading of 
the film on hydrophilic domains. 

To analyze the effect of the length of the stripes, lh, we 
have carried out additional simulations at fixed / while 
varying lh and have measured the corresponding growth 
rate. We plot these results in Fig. SI where we observe 
that the growth rate decreases systematically as lh is in- 
creased. Accordingly, the critical fraction f c {lh) is a de- 
creasing function of lh- 



Analytical model for growth on transverse stripe patterns 

In order to gain insight into the growth of the contact 
line on transverse patterns, we will analyze the inter- 
play between hydrophilic and hydrophobic stripes using 
a kinematical argument. 

Let us denote Ati and Ait the time intervals in which 
the leading and trailing edges sweep a distance Al = 
l<t> + lh, respectively. For a mean steady growth of the 
interface, each of these intervals can be decomposed as 
the sum of the time that the interface spends on a hy- 
drophilic and a hydrophobic stripe. We thus write 



At, 



Aif 



At* 



and 



At t = AiJ 1 + Att > (14) 



where superscripts h and <fi correspond to the hydrophilic 
and hydrophobic domains, respectively. 

For sufficiently long hydrophobic stripes, the thin film 
should grow as an unperturbed finger, with leading and 
trailing edge velocities vi and Vt, respectively. Disregard- 
ing the transient relaxation to these values, it is reason- 



able to assume that the contact line changes its velocity 
instantaneously to vi and v t as soon as it comes into con- 
tact with a hydrophobic stripe. From this assumption, it 
follows that the residence time of the leading and trailing 
edges of the film in a hydrophobic stripe are 



Atl 



l,b 



and 



Att 



'</> 



For the hydrophilic stripes, we can write similar ex- 
pressions for the residence times of the leading and trail- 
ing edges of the interface in a hydrophilic stripe, 



Atl = k 



and 



The velocity v\ corresponds to the mean velocity of the 
leading edge on a hydrophilic stripe, while is defined 
as the velocity of the jump that the trailing interface 
undergoes on a prewet hydrophilic stripe. 

To estimate , we first approximate the volumetric 
flux supplied by the finger to the hydrophilic stripe as 
qi ~ hfXWvi, where hf and \W are the thickness and 
width of the finger in the hydrophobic stripe, respectively 
(see Fig. [5]). The flow rate qi sustains the spreading of 
the leading edge, which covers a volume Vi ~ hflhiihW 
of the hydrophilic stripe before it moves over the next hy- 
drophobic stripe. The volume Vi depends on Wh, which 
is the lateral fraction of the stripe that the film covers be- 
fore the leading edge advances over the next hydrophobic 
domain. In writing the expression for Vi we have assumed 
that the film thickness does not vary significantly during 
spreading. The residence time of the spreading process 
is therefore Ati h = Vijqi, which gives a mean spreading 
velocity 



h 



Vl 



At, h 



w h 



-vi- 



(15) 



We next estimate the jumping velocity of the trailing 
edge, Vf. This velocity corresponds to the ratio between 
the length of the hydrophilic stripe, lh, and the waiting 
time to observe the jump in the contact line position, 
Att ■ Let ht be the thickness of the trailing edge of 
the film when it touches the hydrophilic stripe, which is 
already wet and has a local thickness h = ht° , as shown 
in Fig. [2] After contact, the cross section of the film in 
the x — z plane will evolve as A t h ~ hlh = h t °lh + qtt, 
where the flow rate per unit length in the y direction, 
q t , obeys q t ~ v t h t . We expect that A t h grows until the 
thickness of the film is h = h t , from which it follows that 
the corresponding waiting time is 



AU h = 



(h t - h t °)h 
v t h t 



(16) 



We thus obtain an estimate for the jumping velocity of 
the trailing edge, 



(17) 



8 



where h t = h t Q /h t . 

We now return to Eq. (|14[) , which we write in terms of 
the mean velocities of the leading and trailing edges, vi 
and Vt, i.e., 



l<j> + lh 
vi 



1 JL 

h 



I'l 



Up 
vi 



and 



l<b + lh h 



Vi 



v t ' 



l<t> 
v t 



(18) 



From Eq. (|18p we obtain the mean velocities of the inter- 
face as a function of the local velocities and the fraction 
of hydrophilic stripes, 



vi = 



1 + / 



and 



vt 



l + .f 



>-t h 



(19) 

where vi and vt can be measured from the single finger 
solution, while u ; and v^ follow from Eqs. (fT5")) and (|T71) . 
For the leading edge velocity we obtain 



vi 



vi 



1 



/ ^ 



1) 



(20) 



which is proportional to the leading edge velocity of the 
finger and depends on the fraction of hydrophilic stripes, 
as well as on the degree of transverse spreading, which is 
contained in the dependence on Wh- 
For the trailing edge velocity we find 



vt = 



Vi 



1-htf 



(21) 



which is controlled by the local thickness of the film on 
the prewet hydrophilic stripe. 

Equations (f2T)|) and (|2"Tj) can be used to obtain the 
growth rate of the contact line, 



vi - v t 



L'l 



Vi 



1) 1 - fhf 



(22) 



In order to verify the validity of Eq. (f22j) . values of A, Vf, 
Vt, Wh, an d h t have to be provided. The finger width, 
A, and the finger velocities of the leading and trailing 
edges, Vf, and Vt, are known from the finger growth on a 
homogeneous hydrophobic substrate. For our simulation 
parameters, vi = 1.1, Vt = 0.6, and A ~ 0.3 as measured 
from the finger solution (/ = 0). On the other hand, ht 
and Wh are geometrical parameters that can be measured 
from simulations. We find that h t does not depend on 
the length or on the fraction of hydrophilic domains and 
we measure ht = 0.33. In contrast, Wh shows a depen- 
dence on lh, which is essentially controlled by the compe- 
tition between the forward motion of the film the lateral 
spreading caused by the hydrophobic covering. For long 
hydrophilic domains, the film spreads to a larger extent 
before the leading edge arrives to the next hydrophobic 
domain. Therefore, Wh increases with lh- This introduces 
the dependence of the growth rate on lh ■ We do not have 
a prediction for Wh as a function of lh- Instead, we mea- 
sure these values from our simulations. For the runs per- 
formed we find Wh = 0.58 for lh/x c = 1.3, Wh — 0.60 for 
lh/xc — 4.0, and Wh = 0.64 for lh/x c = 6.7. 



Lattice-Boltzmann 
Theoretical prediction 




w h 



FIG. 5: Critical fraction, f c , as a function of the fraction of the 
hydrophilic stripe covered by the leading edge, Wh- The frac- 
tion Wh increases with the length of the hydrophilic stripes, 
lh- Lattice-Boltzmann simulations (solid symbols) show that 
the critical fraction decreases substantially with increasing 
Wh, i.e., with increasing lh- This behavior is captured by a 
kinematical model for large Wh, depicted by the solid line (see 
text). 



In FigureHJwe show a comparison of Eq. (|2"2")l with sim- 
ulation results. As depicted in the figure, the theoretical 
prediction captures the main dependence of the growth 
rate of the contact line, L, both on / and, through Wh, on 
lh- For a fixed hydrophilic domain length, lh, L decreases 
with increasing / due both to a decrease of the leading 
edge velocity, as shown by Eq. ([20]) . and by an increase 
of the trailing edge velocity, as predicted by Eq. (f2~Tj) . On 
the other hand, for fixed /, increasing lh has the effect 
of decreasing the growth rate. This effect can be traced 
back to Eq. (jT5"j) , which shows that increasing Wh has the 
effect of slowing down the leading edge on a hydrophilic 
domain. 

The theoretical result given by Eq. (|2"2")l predicts a crit- 
ical fraction of hydrophilic stripes, f c < 1, above which 
the growth rate, L, vanishes. Such a saturation emerges 
both from the decrease of the leading edge velocity and 
the increase of the trailing edge velocity, an effect that 
is intimately linked to the fraction of hydrophilic stripes 
imposed to the substrate. The critical fraction follows 
from setting L = in Eq. (|22|) . and reads 



fc = (vi ~ v t ) 



vih t 



Vi 



w h 
A 



(23) 



Here, the relevant dependence is f c (wh), given that Wh 
can be controlled by choosing the length of the hy- 
drophilic domains lh- 

In Fig. [S]we show a plot of f c as a function of lh and Wh- 
We have included f c data for four additional values of lh', 
l h /x c = 0.4, l h /x c = 0.7, l h /x c = 10.1 and l h = 13.4. 
Corresponding values of Wh arc Wh — 0.33, Wh — 0.45, 
Wh = 0.80 and Wh — 0.95, respectively. We also show the 
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theoretical prediction given by Eq. (|23[) . which accurately 
captures the decrease of f c as Wh increases for large Wh ■ 
An interesting limit corresponds to short hydrophilic 
domains, which we can include in our theory. As shown in 
the figure, the critical fraction reaches a saturation value, 
which originates from a weaker lateral spreading. In this 
limit the lateral size of the film in a hydrophilic domain 
remains close to the width of the finger, i.e., Wh ~ A. 
Still, the hydrophilic interaction causes a decrease in the 
local contact angle. As the leading edge encounters the 
next hydrophobic boundary, the capillary ridge grows, 
until the contact angle reaches the advancing contact an- 
gle of the hydrophobic stripe. The relevant timescale is 
given by the filling of the capillary ridge, which we ap- 
proximate as a spherical cap of radius R ~ hf/2 and vol- 
ume Vi = 2ir(hf /2) 3 /3. The volumetric flow coming into 
the sphere still obeys qi = hfXWvi. Adding this contri- 
bution to the residence time, we find the sprcading-filling 
velocity, 



vf(w h -> A) 



6nXWl h T A 



(24) 



Eq. (f2"4"]l can be used to obtain the critical fraction in 
this limit, which reads 



fciw h ->• A) = (vi - v t ) 



vih t 



Vt 



GnXWlt 



(25) 

where Zj£ is the crossover hydrophilic stripe length below 
which filling dominates over spreading. The dotted line in 
Fig. [5] shows the value predicted by Eq. (|25|) for l\jx c ~ 
0.4, which is a typical stripe length for which no lateral 
spreading is observed. As depicted in the figure, this 
limiting value reasonably captures the saturation of the 
critical fraction for small w^. It is important to remark 
that the estimate of the leading edge velocity given by 
Eq. (J2U is expected to be valid only as Wh — > A; for 
larger Wh, the shape of the leading ridge is expected to 
depend on Wh- Still, the filling timescale is subdominant, 
as shown by the good agreement of simulation data with 
Eq. (j2"3")l , and we thus disregard it. 

The opposite case, corresponding to Wh 1, is an 
interesting regime, as in this case the imposed pattern 
has the smallest critical fraction of hydrophilic stripes 
for which saturation is expected. In this limit the leading 
edge spreads completely in the transverse direction on a 
hydrophilic stripe, a situation favored by long hydrophilic 
domains. In Eq. (f2"3"|. this is equivalent to to setting 
Wh = 1, from which it follows that the minimum possible 
critical fraction is 



/ r °° = (vi - Vt) 



vih t 



1 v t 



(26) 



For the simulation values used in this work we find 
f£° ~ 0.3. In Fig.(5Jb) we plot this value, which is a good 
estimate of the limiting critical fraction as suggested by 



Lattice-Boltzmann results. While both simulation and 
theory predict a slow final decay towards the limiting 
critical fraction, there is a rapid variation for intermedi- 
ate stripe lengths, as shown in Fig. [SJ For example, for 
h/xc — 13.4 we obtain f c ~ 0.4. This value of the stripe 
length should be compared to the the typical finger spac- 
ing, which for our simulations is A max /x c = 8. There- 
fore, one expects critical fractions significantly smaller 
than unity for stripes whose length is comparable to the 
lengthscale of the unstable film. 



V. CONCLUSIONS 

By means of Lattice-Boltzmann simulations and ana- 
lytic approximations, we have studied the effect of hetero- 
geneous hydrophilic-hydrophobic substrates composed of 
transverse striped patterns on the dynamics of forced thin 
liquid films, and have demonstrated that the film growth 
can be controlled by the interaction of the film with a 
prescribed hydrophilic-hydrophobic pattern. 

We have focused on a scenario where the unstable 
contact line gives rise to sawtooth and finger structures 
on homogeneous hydrophilic and hydrophobic substrates, 
respectively, and have examined how the growth of the 
contact line is altered by the chemical pattern. To 
this end, we have considered patterns where the typi- 
cal lengthscale is comparable to the most unstable wave- 
length of the contact line, A max . 

For the transverse striped patterns considered, we have 
shown that above a critical fraction of hydrophilic stripes, 
f c , the saturation of the film is achieved. The mechanism 
that causes such a saturation originates from transverse 
spreading of the leading edge of the contact line on the 
hydrophilic domains and from the growth of the film at 
the hydrophilic-hydrophobic boundary so as to reach the 
advancing contact angle. These processes contribute to 
a reduction of the velocity of the leading edge. As this 
occurs, a film of fluid is left on the hydrophilic domains, 
thus increasing the velocity of the trailing edge as soon 
as it comes into contact with a hydrophilic stripe, which 
already has been wet. The trailing edge jumps every 
time it finds a wet domain, with the corresponding gain 
in the trailing edge velocity. The global outcome is that 
for sufficiently large fractions of hydrophilic stripes, the 
net growth rate of contact line vanishes, independently 
of the details of the dynamics of the contact line. 

We have captured this effect with a kinematical model, 
assuming local relaxation of the leading and trailing edge 
velocities to 'hydrophilic' and 'hydrophobic' values. Hy- 
drophobic values have been approximated as the veloci- 
ties of the leading and trailing edges of an unperturbed 
finger, while the hydrophilic velocity values have been 
estimated by mass conservation. The global growth rate 
of the contact line, L, has been determined by the dif- 
ference of the mean velocities of the leading and trailing 
edges of the contact line. By estimating these velocities 
on kinematical grounds, we have obtained a prediction 
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for the global growth rate, which is given by Eq. ([22]). 
We have estimated the critical fraction, f c , at which sat- 
uration occurs indirectly as a function of the length of 
the hydrophilic domains, given by Eq. (|23[) . For long hy- 
drophilic domains, in which transverse spreading is more 
efficient, we find a good agreement between the simula- 
tion results and the theoretical prediction. In this limit, 
we have estimated the minimum fraction of hydrophilic 
stripes for which saturation is expected, which for the 
simulation parameters considered in this work is as low 
as f c ~ 0.3. However, experimentally feasible domain 
lengths, comparable to the length scale of the growing fin- 
gers, give critical fractions of about f c = 0.6 to f c = 0.7, 
still considerably smaller than unity. 

Our theoretical prediction, given by Eq. (|2U)) can 
be used to estimate the critical fraction of hydrophilic 
stripes above which the saturation of the contact line is 
expected. To examine the reliability of our approach in 
a real system we consider a potential microfluidic device 
where a water film of thickness h c = 10 /im is forced 
on a micropattcrn in the presence of air. The air-water 
pair closely reproduces the viscosity contrast used in the 
simulations. Our simulations show that the instability 
is effectively suppressed using a sharp contrast in the 
wetting properties between the domains forming the mi- 
cropattern. For example, a front moving over a Poly- 
methyl methacrylate (PMMA) substrate would have a 
relatively large contact angle (9 e « 74°), where fingers 
are expected to develop Q . In that case, the low- wetting 
angle domains could be formed using glass [6 e « 25°) or 
a gold coating (0 e ss 0°) [Hj]. 

Microfluidic devices operate over a wide range of veloc- 
ities, covering from /im/s to cm/s [33j . Here we consider 
an injection velocity of u = 1 mm/s as a representative 
example. Under these circumstances, the required body 
force to drive the film follows from the force balance 
far from the contact line, f x = 3r]iu/h^. ~ 10 4 Pa/m. 
Such body forces are of the order of the gravitational 
force density for water and can be accomplished by sim- 
ple column devices and accurately controlled by syringe 
pumps. From our simulation measurements we estimate 
the leading and trailing edge velocities of steady finger 



propagating at the imposed speed on a hydrophobic sub- 
strate, vi « 1.1 mm/s and vi ~ 0.6 mm/s. Our prediction 
shows that the film will saturate at a different fraction 
of hydrophilic coverage depending on the width of the 
hydrophilic stripes, lh, which determines the fraction of 
lateral spreading, wt ■ A sensible choice is to set lh to the 
order of magnitude of the lengthscale of the instability, 
%c = /i c (3Ca) -1 / 3 . Thus, for lh — x c ~ 10 2 /im, we get 
the spreading fraction Wh ~ 0.58. Using this value of 
Wh and comparing with Fig. 4, the film is expected to 
saturate at a critical fraction f c « 0.6. 

The previous example shows the experimental feasibil- 
ity of the mechanisms proposed in this work for a reg- 
ular liquid under typical conditions found in microfiu- 
idics. However, the fact that the relevant lengths scale 
with the film thickness (expressed here by the choice of 
units to report our results) indicates that a proper up- or 
down-scaling of the system size and the applied forcing 
can provide a substantial flexibility to induce analogous 
phenomena in different liquid/gas or liquid/liquid mix- 
tures at scales below the capillary length. This opens 
the possibility of a better control of the fingering instabil- 
ity in controlled environments by substrate heterogene- 
ity. Our approach relies on chemical patterning, which 
has already been explored experimentally, even at the 
microscale. We therefore expect our results to motivate 
further experimental investigations, with potential im- 
pact in microfluidic technologies. 
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